Extracting Quasiparticle Lifetimes from STM experiments 



O 



5-H 

Oh- 



o 



o 



Sumiran Pujari 

Department of Physics, Cornell University, Ithaca, New York 14853-250fc\ 

Based on Quasiparticle interference(QPI) around a point impurity, we demonstrate an analysis scheme that 
extracts the lifetime of a quasiparticle by using the local density of states(LDOS) data around the impurity in 
a Scanning Tunneling Microscopy(STM) experiment. This data analysis scheme would augment the Fourier- 
Transform Scanning Tunneling Spectroscopic methods which provides us with the quasiparticle dispersion. 
Thus, point impurities can be used as probes to extract quasiparticle lifetimes from STM experiments and this 
would complement other experimental methods such as Angle Resolved Photo-emission Spectrocopy(ARPES). 
We detail how the scheme would apply to metals and superconductors. 



Scanning Tunneling Microscopy(STM) has revolutionized condensed matter research by providing us with unprecedented 
detail on the local real space electronic properties of the sample under investigation. But even more remarkably, it has been 
shown that momentum space properties of the sample can be extracted through the application of Fourier Transform Scanning 
Tunneling Spectroscopy (FT-STS) III]. Using FT-STS one figures out the dispersion of the underlying well-defined quasiparticles 
or carriers. 

In this paper, we aim to extend the domain of momentum space properties that can be extracted using STM. The central 
result is the demonstration of a data analysis scheme that would give us the lifetimes of the charge carriers in a sample as a 
function of momentum(and energy) from data collected in an STM experiment. Previously, ARPES is the tool that has been 
used successfully to extract lifetimes of carriers in a sample by measuring the one-particle electron spectral function directly 
in momentum space. Extracting lifetime information from STM - a real space probe - thus would add value by providing an 
independent method that complements and checks the ARPES method. Previous attempts at reconciling lifetime broadening 
effects on STM data mainly consist of writing down viable fitting forms for the lifetime function that fit with the STM data, 

rather than extracting it out of the data directly like one does in an ARPES experiment by quantifying the width of the peaks in 

r~i 

ARPES spectra (See e.g. [2]). In the context of metals/Fermi liquids, Ref. 130 have fitted STM data on Silver and Copper with 
a model for thermal broadening of the electrons B4|]. Refs. Q5J] and ||6|] are prominent examples in the STM phenomenology of 
high temperature superconductors. 

We start by describing the scheme in the simpler case of normal metals. We imagine the system to be composed of Landau 
quasiparticles described by a propagator of the form 



G (k;uj) = , 1 — (1) 

^ ' uo — ir]{k,uo) — eik) 

Q\ ; 

Self-energy processes - e.g. due to electron-electron interaction as in a Fermi Liquid or through scattering off a bosonic mode like 
phonons - lead to a finite lifetime for the quasiparticle and this is formally taken care by the imaginary term in the denominator 
of Eq. (Q]), rj{k, uo). Also, the real part of the self-energy shifts the chemical potential and we assume that the dispersion term 
e{k) is this shifted dispersion [7]. Our aim is to extract r](k, uo) from STM data. We assume the knowledge of the dispersion 
e(k) either via FT-STS on the same data or through an ARPES experiment. 

Apart from the quasiparticles, let us imagine there to be a point impurity in the system, say at origin, which scatters the 
quasiparticles. In a real situation, we are imagining there to be a dilute amount of impurities in the sample so that multiple 
impurity scattering is not important. The impurity problem is solved via the T-matrix approach [8], and the real space impurity 
^ | scattered electron propagator is given by 
■ 

_ G(r, t'\uj) = G (r, r'; uj) + G (r, r imp ;uj) • T(uj) • G (r imp , r'; uj) (2) 

where Go(r, r uj) = (2tt/L) 2 Go(k; uj)e lk ^ r ~ r ) = Go(R = r — r'\ uj) is the free electron propagator and the impurity 
effect is captured by the so-called "T-matrix" T(uj), which is given by T(uo) = U/(l — U Go(ri mp , ri mp ; uo)) where U is the 
impurity strength. It is in the second term of the above equation that we have QPI which is utilized in FT-STS. 

We will quickly review the key notions underlying FT-STS, since our method also utilizes QPI. STM measures LDOS as 
a spatial map over the surface for a range of energies. The LDOS n(r;uo) is proportional to imaginary part of the real space 
propagator, i.e. 

n(r;uo) = Im[G(r,r;uo)]. (3) 
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FT-STS 's main operating principle is that the peaks in the Fourier transform of LDOS map at a particular energy are at wave- 
vectors which connect pairs of points on the e(/c)'s contour at that particular energy for which the joint density of states is 
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FIG. 1 : In this figure, we demonstrate the effect of Kramers-Kroning to an example LDOS where we limit the integral by a finite cut-off, 

Re[G(r,r;uj)] = P f_ 

a (J-x) - ^ ne example LDOS (see inset) is for a nearest-neighbour hopping model at half-filling, Go(k;uj) — 
(uj — i0.lt + 2t(Cos[k x ] + Cos[/c y ])) _1 and t — 1. Around the Fermi energy (cj = in this case), we see that even for |A| = 3, Re[G] agrees 
well upto around \uj\ = 1. One can quantitatively show that this error is at most Log\ | « 2\u/A\ in units of n and we do much better 
than that(see Supplementary). 

maximum. This can be understood by looking at the Fourier transform of the interference term in Eq. (|2) (see Eq. (1) of 
and the following paragraph). If the quasiparticles have finite lifetimes, the resultant effect in FT-STS will be a broadening of 
the FT-STS peaks(which are seen in experiments, e.g. [10]). Moreover, the "shapes" of these FT-STS peaks contain information 
about the momentum dependence of the lifetime r](k; uj). It seems that extracting the /c-dependence of r](k; uj) from the FT-STS 
method is a hard task because, apart from other possible broadening factors like inhomogeneity (e.g. STM on cuprates), one has 
the difficulty of deconvolving the output of FT-STS - the QPI term is a product in real space - without the prior knowledge of 
rj(k; uj). Instead we will work in real space, our main tactic being to extract Go(R; uj) from QPI, and STM data is most suited 
for this. 

We now list down the main steps of the analysis scheme and in what follows we give their essential technicalities along with 
pictorial demonstrations. In the Appendix, we include further technical details and proofs required in those steps. 1) From 
LDOS/n(r; uj) map, we construct a G(r, r; uj) map. 2) Once we have the G(r, r; uj) data, we "invert" Eq. © in order to extract 
Go(R; uj)- To invert Eq. ©, we need 2a) a way to find Go(R = 0; uj) and 2b) a way to find the correct phases of Go(R; uj). 
Once this is done, we Fourier transform to get Go(k;uj) and, thence, r](k;uj). We show results of this method for various 
cases of dispersion and lifetimes. Then, we discuss what kind of data sets are desirable and how the method extends to the 
superconducting case. 

The first step of the analysis method is to convert the LDOS data to G(r, r; uj). This will be achieved through a Kramers- 
Kronig relation the propagator satisfies, Re[G(r, r; uj)] = P J n(r, r; x)/(uj — x) where the principle value integral is over the 
real line. Since the LDOS is nonzero only within a finite bandwidth [11], this integral is over a finite range of energies. In 
general, in a real experiment one might have information only over part of the bandwidth in which case, we can definitively 
apply this method only to an energy range that is well within the measured energy range, where even the incomplete spectrum 
can be fruitfully used as demonstrated in Fig. [T] This is very often the/one of the interesting energy ranges(e.g. around the Fermi 
energy for metals or the nodal energy for cuprates). We can also apply some form of extrapolation to construct LDOS data over 
the full bandwidth 11211 . Kramers-Kronig has been applied successfully to other spectroscopies, e.g. Electron Microscopy (see 
fl3ll ). thus giving us reason that they be applied to STM data as well. 

We now discuss the second step : how to invert Eq. © at a fixed energy. We are only concerned with r — r' '. We set 
Timp — 0. The first step is to find out the first term on the right hand side of Eq. ©, Go(R = 0;uj). This will be done through 
a minimization procedure where a cost function would penalize incorrect guesses for Go(0; uj). Given a Go(0; uj) guess(which 
is independent of R if the free propagator is that of a translationally invariant system), we can solve for T(uj) by solving Eq. 
for r = r' = ri mp (Furthermore, we can calculate the impurity strength U from T(uj)). Once T(uj) is known, we can solve for 
Gq(R; uj) as 
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FIG. 2: Above is shown \Go (R; uj) | on a 200-site window around an impurity extracted with various start guesses for Go (0; uj). 

In Green's function theory, one can show that the magnitude \Gq(R\uj)\ monotonically decays to zero for large R (expo- 
nentially in R in one dimensions and as square root of R in two dimensions, see Supplementary) for dispersion that have 
convex energy contours. We demonstrate this effect in ID and also show the effect of incorrect Go(0;cj) on extracted 
\Gq(R;uj)\ in Fig. [2] We see how an incorrect guess for Go(0;u) spoils the monotonic decay of \Go(R;uj)\. The reason 
for the deviation from monotonicity is as follows : Given our (incorrect) guess of |Go(0; uj)\, we can decompose the incorrect 
G (R;uj) as G c orrect (R;uj) + G e rror where G e rror is a constant. Therefore, \G (R;oj)\ = \G c orrect (R;uj)\ + \G e rror \ + 
2\Gl orrec \R- uj)\\Gl rror \ x Cos(Arg[G^ orrect (R; uj)} - Arg[G^ rror }), and it is the final cosine term in the above expression 
which spoils the monotonicity even for large R. Moreover, the |G?g rror | term would also not let the Green's function decay to 
zero as r — » oo. This motivates a minimization using a cost function that penalizes deviation from the smooth decay of extracted 
\Go(R; uj) I for finding the correct |Go(0; uj) | J3. A good start guess for Go(0; uj) is to take a spatial average of G(r, r; uj) over 
the whole data set around the impurity. One can show that the error in the guessed Go(0; uj) is 1/L 2 (l/L d in d dimensions, see 
Supplementary) suppressed compared to the guessed Go(0; uj), and if the window were infinite, the spatial average of G(r, r; uj) 
would exactly equal Go(0; uj). 

With the correct Go(0; uj), we still get Go(R; uj) only up to a phase of tt. Capturing this phase is crucial to get the correct 
Go(k; uj) upon Fourier transforming. To get the correct phase, we start with the observation that the phases have to be smooth 
and well-behaved as a function of R because Go(R;uj) is differentiable with respect to R lfl5ll . We use this property to fix the 
phase of the square of Go(R; uj), i.e. we select that branch of the argument function when evaluating the phase of Go(R; uj) 2 
which maintain the aforesaid smoothness. We start by making a spatial list of the phases as given by the Arg(z) function which 
restricts the phase obtained to one branch of the Argument function. Then, we start at R = 0. As we move away from the 
origin, we multiply phase factors of e l2m7V to Go(R;uj) 2 = \Go(R'- ) uj) 2 \e %( f )principal for all R, the m's being so chosen that if 
\R'\ > \R\ then (j)' principa i + 2irm / > (^principal + 2irm. Once that is done, the phase of Go(R; uj) is just half that of Go(R; uj) 2 . 
We demonstrate the working of this phase reconstruction method in Fig. [3] (see Supplementary for a flowchart of the method). 

With the correct phases, we are now ready to Fourier transform the extracted Go(R; uj) to get Go(k; uj) and r](k; uj) with 
our knowledge of e(k). Moreover, the extracted Go(k;uj) also has to satisfy the exclusive momentum dependence of uj — 
Re[Go(k; uj) -1 ]. In Fig. IH we show how this method performs with and without error and we see that it performs well for error 
magnitudes less than 0.25 %. For the panels Fig. |4]a-d, the form of r\ had no momentum dependence, and this kind of fitting form 
has been proposed in [6] for Cuprates and has been theoretically discussed in 111 611 . In general, we expect the lifetime function to 
have few (low) harmonics of k similar to the dispersion. Thus, our analysis scheme would serve the purpose of finding the most 
general r](k;uj) that is consistent with STM data. We can extract an approximate analytic form for r] by doing a least-squares fit 
of the extracted r\ to a function of k containing a few harmonics in the Brillouin zone. The approximate analytic form can then 
be compared to theoretical proposals. 

At this point, we comment on what kind of data sets would be ideal for such an analysis. In Fig. [51 we show an example 
of data set seen in a real experiment. We show how it is similar to a theoretical data set(calculated numerically) which has a 
lifetime broadening. Thus, we would expect that if we observe a few of the "Friedel oscillation"-like rings around the point 
impurity, this analysis scheme should work. Moreover, if FT-STS applied to a single point impurity data shows reliable QPI 
peaks, then we believe that the data set would have good enough spatial resolution to resolve the momentum dependence of 
lifetime r\ to the same momentum resolution as that of the FT-STS results. We can improve on this by taking an average over 



FIG. 3: In the panels above we show the phases of Go(R; uj) (a) and b)) and Go(R; uj) 2 (c) and d)) in the first quadrant of size 30x30 lattice 
sites around an impurity at a fixed u(= — £). a) and c) show the phases as evaluated by the Arg(z) function restricted to one branch, b) and 
d) show the smooth phases as reconstructed using the reconstruction algorithm. The ratio of phases in b) and d) is identically two over the 
whole quadrant, even though the ratio of phases in a) and c) does not behave in such a regular manner. e(k) = — 2t(Cos[k x ] + Cos[k y ]) and 
rj(k, uj = —t) = O.lt in this example. 

data sets around multiple point impurities to improve signal to noise. Now, we will sketch how this method of analysis can 
be extended to superconducting case using d-wave superconductors (pertinent to Cuprate phenomenology) as our example. In 
Nambu's two component notation, the free superconducting propagator looks like 



where e(k) is the bare dispersion and A(fc) is the (d-wave) gap of the cuprate in question. These are assumed to be 
known(through other experiments). As before, we want to determine the electron/hole lifetime. The first simplification is the 
relation rjh(k,u) = —r] e (k,—oj). The proof of this relation in outlined in the supplementary information to this manuscript 
and it follows by showing E22(— w — id) = —En (a; + id). This relation implies G$(R, 00)22 = —Go(R j —oj)u and 
Go(R, o;)i2 = Gq(R, —uj)i2- Now, as before, we imagine there is a point impurity which result in a two-component T-matrix. 
One can show that this T-matrix has no off-diagonal entries(for an ordinary potential impurity) since Go(R = 0,gj)i2 = 
owing to the d-wave symmetry of the gap function. One can further show that T22(— uj) — —Th(oj) and, resultantly, 
G(r,r;uj) 2 2 = -G(r,r; -uj)u. For r = r imp , we have G n = G ,n + ^11^0,11 and G22 = G 0y 22 + ^22Gg j22 . Using 
G(r, r; , 00)22 — — G(r, r; , — cj)h, we can thus determine Tn and T22 given a guess for Go(R = 0, oo)u (which will again be 
determined by demanding the monotonicity of Gq(R, )• For r =^ ri mp , we have 



G (k;oj) 



-1 



00 — irj e (k, oj) — e(fc) 
A(fe)* 




(5) 



Gn(r,r;o;) = G (0, u)u + Tn(cj)G (r - r irnp , oj)^ 

+T 22 (^)G (r - r imp ,oj)i2Go(r - r imp ,cj) 2 i 

G 2 2(r,r;oj) = G (0, oj) 22 + T 22 (^)G (r - r irnp: oj) j 2 

+Tn(o;)Go(r - r imp ,cj)i 2 G (r - r imp ,cj) 2 i. 



(6) 



(7) 



Again using G(r, o;)22 = — G(r, — oj)h, now with the knowledge of Tn and T22, we can solve the above equations for Go(R, oj) 
upto a phase of tt (which we reconstruct as before) at each oj for all r in the dataset, following which we Fourier transform to 
extract Go(fc, oj and rj e (k, oj). 




FIG. 4: In these figures we are plotting |(cj + irj(k,uj) — e(k))~ 1 \ as a function of k over a Brillouin Zone (0, 2tt) x (0, 2tv) at a fixed 
uj(= —t for the above plots). In a) we show the input form resulting out of our choice of input for 77, where e(k) = — 2t(Cos[k x ] + Cos[k y ]) 
(nearest-neighbour hopping) and 77(A), uj = —t) = O.lt; in b) we show the form extracted using the proposed analysis scheme when no noise 
was added to the STM data calculated numerically. One sees the limitation in momentum resolution in the form of "blockiness" introduced 
by having a finite window. This "blocky" momentum resolution gets better or worse with greater or smaller window sizes. In c) and d) we 
show the results of the analysis scheme to data with 1% and 0.05% Gaussian errors added respectively. We have done similar analyses for 
different energy values and different forms of 77 and in e), f) and g), we show the corresponding results for e{k) — — 2t(Cos[k x ] + Cos[k y ]) — 
4(0.2t)(Cos[k x ] * Cos[k y ]) (nearest and second-nearest neighbour hopping) and r](k,uj = —t) = 0.25£ + 0.1t(Cos[k x ] + Cos[k y \) as 
another example. 
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FIG. 5: In this figure we compare an experimental data set on In As surface(taken from H7I1 with due permission from APS and the authors) 
with a numerically calculated LDOS data set with error 0.5 % added. This figure serves to illustrate that there exist data sets, perhaps within 
operable error range, to which the scheme can potentially be applied. 



In conclusion, we demonstrated an analysis scheme which holds promise to extract lifetimes from STM data in various 
systems ranging from metals and semiconductors to strongly correlated compounds to superconductors. Some final remarks are 
in order. We demonstrated the proposed analysis scheme in case of a point impurity, but it can be extended to the case of an 
extended impurity too. The resulting complication will be that now we would have to guess more than just Go(R = 0; uj) (e.g. 
if the impurity extends over two adjoining sites n and r2, then G(r, r; uj) will be a function of Go(0; uj) = Go(r±, ri; uj) — 
G^{t2\T2\uj) and Go(l; uj) = Go(rl, r2; uj)). This scheme is inherently local, where we would be analyzing data around a 
single impurity. Thus, it would really utilize the local information that STM affords us with. There have been other examples 
of data analysis done on STM data previously to extract localinf ormationf 1181. Jl9ll ) . In this sense, we would do better than 
ARPES where the signal is averaged over an area of the sample equal to the beam size, iff the STM experiment has good signal 
to noise. Similarly, we can overcome inhomogeneity issues for dirty systems, in which case we would concentrate this analysis 
on a homogeneous patch similar in spirit to Hudson's analysis lfl9ll and to a previous work [I20I1 . 
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suggestions of C. L. Henley, especially regarding the checks on robustness of the method to error and phase fixing, and a critical 
reading of the manuscript. I also thank him and H. J. Changlani, Milan Allan and J. C. Davis for helpful discussions. 



APPENDIX : SUPPLEMENTARY INFORMATION 
Limit on the Error introduced by Kramers-Kroning Integration 

The Kramers-Kronig relation relates the real part of a Green's function to the imaginary part as follows 

^ r^/ \ -1 1 „ f°° 7 Im\G(r,r:x)] „ f°° 7 n(r:x) 

Re\G(r,r;uj)} = P dx ) { = P dx . v 7 \ . (8) 

7T J_ 00 [uj - x) J_ 00 (uj - x) 

If we limit the integral by cut-offs A + and A_ , then the error E introduced is 



J-00 - X) Ja + 



(uj — x) 



Even if we were to make the really bad approximation that n(r;x) = n(r;uj) for all x (and this is a really bad approximation, 
since n(r; x) decays to zero as x — » 00), we get 

E <n(r-uj)Log ^~ K -\ (10) 
\UJ-A + \ 

since n(r; x) > for all x. Thus, if cj/A_ << 1 and uj / k+ « 1, then the error in Re[G(r, r; uj)] is less than n(r, uj) * + 
-ft?- + Logj^). Since, Re[G(r,r;uj)] and n(r;uj) carry the same dimensions and n(r;uj) < 1, this is at most an 0(|^|) error. 
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Proof of Monotonic decay of Go (R; uj) for large R 

In two dimensions, Go (R) on a lattice is given by the formula 

G (r,r irnp ;uj) G (f - r imp ;uj) = G (R;uj) 



lim \ / dk 

Je 



ik.R 



5^0+ (2tt) 2 J BmZm uj + i5- e(k) 



— 7 / dk x e ikxRx / dfc, (11) 

^) 2 J-n J-* V u + i5-e(k) 



(2k)*J_ v J_ n y u + i5-e(k) 

Let us look at the k y integral for a particular k x . The denominator vanishes for certain values of k y thus motivating the conversion 
of the k y integral to a contour integral. The mapping z = e %ky achieves the conversion which also maps the integral from — ir to 
7r to a contour integral over the unit circle. The periodicity of the integrand over the zone ensures the analyticity of the resulting 
complex integrand. Thus, 

~* If* k R I dz Z Ry 

G ^ R)UJ ^ (2tt) 2 7_ 7r ^ e J uc iz uj + i5 - e(k x ,z) 

For a particular uj energy contour and k x , we get two poles (See Fig.). Expanding the denominator around the poles gives us 
u + i6-e(k x ,z) = - hVgy ^ kx \ z- z p (l - j—^j^j ))' The poles z p s are denned by e(k x ,z p ) = uj and hv 9y (uj,k x ) = 
is the group velocity along y direction. We need only worry about the (z — z p ) term in the expansion of the denominator since 
other expansion terms will yield zero residues. From the (1 — ^— — ^ k ^ ) factor in the expansion, we realize that the pole where 

the sign of the Hv 9y is same as the positive S will be "pulled" inside the unit circle while the other pole will be "pushed" out of the 
unit circle. Thus, when we do k x integral, only one half of the uj energy contour (not to be confused with the complex contour; 
to distinguish we will call uj contours as energy contours) will contribute to the integral. In the process, we have converted the 
2D integral over the zone into an integral over part of the energy contour. Filling in the steps, 

dz z Ry 



Go(R;uj) = — i-r / dk x e ikxRx I , , . , 

V J (2^) 2 J_ v fun iz _ ft* fl| >,fc*) 



c lz - g y 1 , 3:; f2 - z (1 2 ^ 



i r 



dk T e ikxRx 2m- 



( 27r ) 2 J-7T ^ hVg y (uJ,k x ) 

^ r e ik x R x e iky p (u,k x )R y 

(2tt) 2 J x hv 9y (uj,k x ) 



i r s,ik(s,uj).R 

= : & ds ^ (13) 

i S gn(5)=sgn(v gy ( W ,s)) | Ve(fc(s, Uj))\ 

where the last step was achieved by converting the element dk x to a parameter s characterising the uj energy contour and we 
integrate over that part of the contour where the sign of S is same as v Qy . 

For large R, i.e. far from impurity, we notice that the phase e tk ( s > UJ )- R varies rapidly and thus stationary phase approximation 
can be applied. The phase factor is stationary at points on the energy contour where the group velocity is along the R direction 
since d(k(s,uj).R)/ds = R.dk(s,uj)/ds = only when R is perpendicular to dk(s,uj)/ds and dk(s,uj)/ds, being the tangent 
to the energy contour, is perpendicular to the group velocity. Therefore, 



G (R;u) = A — — = — e ik dom (R, u) .R 

2™ \Ve(k dom (R,uj))\^ \R\\<Pk{s,u)/ds\ dom{nuj) 

(14) 

where kdom(R,u) is the k corresponding to which the group velocity at energy u is along R and, thus, is a function of R (only 
through R) and u. For a convex dispersion function e(k), we will have only one kdom and thus 

|G (i2;w)|oc— L (15) 
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for large R. This proves the monotonic decrease of Gq{R\ uj) when the lifetime is infinitesimal. When we have a finite lifetime 
due to self-energy processes, the propagator in momentum space looks like G$(k\ uj) = (uj + iS — (c(k) + irj(k, where 
the uj might have undergone a chemical potential shift, and the whole algebra in the above will go through similarly and we will 
get 



\Go(R;uj)\(x —=e ive ( fc dom( iUo)i 1 1 (16) 
In one dimension, we only get the monotonic exponential decay for large R. 



IMPLEMENTATION OF COST FUNCTION FOR FINDING G (R 0; uj) 

As mentioned in the main manuscript, the Cost function for a one-dimensional list of values for \Go(R; uj) | (that is extracted 
given a guess Go(R = 0; uj)) was 

n +(( ^\ ^ \z r +i + z r _i — 2z r \ 2 

where the list is {z r } and sum runs over all 3 -tuples. We generalize this to two dimension by evaluating the one-dimensional 
cost using the same formula for all one-dimensional slices of the two-dimensional data set either along x or y direction. We do 
it this way because the two-dimensional data set is symmetric with respect to interchanging x and y when there is no error. In 
the error- full case, we can pre-process the data set to impose the symmetries of the square lattice. Thus the Cost function is 

We show the profile of this Cost function as a function of Gq(R = 0; uj) guesses for the no-error case (which includes the 
numerical error incurred during two-dimensional Numerical Integration in Mathematica) and error- full cases in Fig. [6] We show 
it as a matrix where the center point(6,6) corresponds to the correct Go(R = 0; uj) and the (7,7)-entry corresponds to average 
over the G(r, r; uj) set (see next section). From point to point, we change the guess by Re[Avg(G(r, r; uj)) — Go(R = 0; uj)] 
along ^-direction and Im[Avg(G(R; uj)) — Go(R = 0; uj)] along y direction. We see that in the no-error case, the Cost function 
has minimum at the correct value of Go(R = 0; uj)]. In case of 0.1%, it does well to within Avg(G(r, r; uj)) — Go(R = 0; uj). 
In case of 0.5%, it starts to seriously deviate and the best guess then would be Avg(G(r, r; uj)). 



A Good Guess for Go(r, r; uj) 

Recalling that the T-matrix equation for scattering of point impurity is 

G(r, r'\uj) = G (r, r'; uj) + G (r, r imp ; uj) • T(uj) • Go(r imp , r'; uj), (19) 

when we take the average of Eq. [19] with r = r' over the window, the two terms on the right hand side average to(in two 
dimensions) 

i^G o (0;a;) = (^) 2 ^G (fc;u;) (20) 

r k 

^^G (r,0;o;)T(a;)Go(0,r; W ) = (^) 4 T( w )^G (fc lw ) 2 (21) 

r k 

(22) 

Thus, we see that the second term is 1/ L 2 (l / L d in d dimensions) suppressed compared to the first term, and if the window were 
infinite, the spatial average of G(r, r; uj) would exactly equal Go(0; uj). For a finite but large enough window, it is a good guess 
forG (0;u;). 



9 



a) 



b) 



53 9. 


23 8 


517. 


192 


4 8 0. 


55 8 


4 55. 


35 7 


483. 


146 


526. 


168 


532. 


907 


481. 


57 8 


479. 


351 


479. 


911 


510. 


442 


50 0. 


72 


497. 


93 9 


45 8. 


983 


4 5. 


2 8 


412 . 


83 


422 . 


243 


419. 


665 


402. 


611 


427. 


5 97 


462. 


664 


514. 


9 64 


50 6. 


821 


4 64. 


84 4 


38 0. 


62 8 


33 9. 


7 4 9 


306. 


183 


315. 


05 


314. 


686 


323 . 


685 


3 91. 


342 


4 50. 


97 5 


517. 


7 2 


4 7 0. 


2 9 6 


4 9. 


64 8 


32 0. 


441 


2 5 8. 


852 


183. 


675 


169. 


546 


193. 


139 


2 67. 


715 


358 . 


4 6 9 


4 64. 


87 


5 7. 


9 4 7 


4 65. 


8 8 8 


411. 


16 6 


30 9. 


8 91 


18 7. 


92 5 


82 . 9017 


45.7886 


86. 9743 


18 8. 


98 


3 7. 


511 


3 98. 


2 8 8 


4 60. 


5 9 9 


485. 


666 


410 . 


571 


313. 


792 


184 . 


867 


48 . 9257 


3.5936 


47 . 1426 


162 . 


098 


278 . 


544 


382 . 


236 


432 . 


65 


521. 


963 


460. 


674 


337. 


781 


223. 


156 


88.8926 


43.0211 


92 .7411 


213. 


409 


300 . 


548 


401 . 


974 


472 . 


436 


576. 


2 


449. 


957 


362. 


884 


272. 


033 


198. 


329 


144 . 


986 


193. 


234 


298 . 


064 


389. 


643 


441 . 


698 


491 . 


246 


566. 


162 


460. 


782 


412 . 


75 


357. 


972 


287. 


743 


254. 


063 


279. 


909 


386. 


979 


487 . 


795 


524 . 


563 


523. 


213 


517 . 


34 


472 . 


211 


448 . 


793 


417 . 


734 


381. 


36 


355. 


206 


376. 


472 


413. 


623 


503. 


829 


582 . 


937 


578 . 


792 


487 . 


028 


471 . 


423 


483. 


924 


498. 


604 


475. 


834 


417 . 


123 


426. 


282 


450 . 


414 


513. 


555 


592 . 


555 


617 . 


897 


581. 


8 65 


5 65. 


535 


55 4. 


34 9 


5 61. 


683 


548. 


3 


556. 


993 


598. 


945 


6 61. 


841 


647. 


13 8 


617. 


003 


619. 


641 


5 90. 


513 


5 6 9. 


113 


552 . 


95 


531. 


317 


531 . 


274 


559. 


229 


595. 


314 


651 . 


2 93 


631 . 


7 01 


63 6. 


191 


64 8. 


2 53 


616. 


95 4 


58 9. 


803 


551. 


7 94 


52 4. 


22 4 


516. 


649 


550 . 


219 


592 . 


995 


60 7. 


17 


613 . 


04 5 


7 02. 


882 


6 69. 


7 4 8 


63 9. 


416 


610. 


02 7 


57 8. 


4 6 9 


52 6. 


832 


499. 


443 


512. 


808 


547 . 


834 


530 . 


7 4 9 


5 82. 


62 5 


710. 


50 3 


645 . 


552 


64 5. 


32 5 


64 8. 


301 


60 8. 


452 


5 67. 


35 


505. 


849 


476. 


692 


483. 


789 


4 8 4. 


80 6 


535 . 


6 


5 90. 


92 6 


637 . 


14 8 


664 . 


955 


700 . 


321 


634 . 


724 


582 . 


189 


514 . 


757 


452 . 


253 


438 . 


86 


454 . 


124 


501 . 


073 


556. 


25 


613 . 


204 


689. 


518 


744 . 


853 


620 . 


959 


541 . 


026 


498 . 


365 


458 . 


653 


439. 


113 


453 . 


436 


519. 


962 


559. 


975 


577 . 


771 


665. 


557 


644 . 


157 


602. 


54 


552. 


118 


502. 


903 


477 . 


294 


470 . 


276 


506. 


134 


523. 


295 


550 . 


928 


574 . 


724 


655. 


917 


635. 


459 


610. 


326 


562. 


764 


530. 


082 


518 . 


215 


522. 


233 


561 . 


17 


569. 


449 


606. 


031 


598 . 


648 


623. 


632 


617. 


58 


618. 


168 


610. 


082 


559. 


48 


538. 


606 


550 . 


546 


586. 


173 


612 . 


755 


604 . 


402 


603. 


033 


610. 


287 


630. 


03 


657 . 


41 


642. 


103 


599. 


747 


575. 


798 


594 . 


722 


622 . 


91 


662 . 


804 


625. 


179 


609. 


922 


950. 


21 


971. 


854 


1009 


.56 


1048 


.68 


1079 


.02 


1104 


.7 


1133 


.77 


1177 


.03 


1186 


.42 


1136 


.86 


1097 


.61 


974. 


27 


993 . 


663 


102 8 


. 7 6 


10 90 


.11 


1111 


.41 


1112 


.89 


1127 


.65 


1174 


.47 


1225 


. 3 8 


113 8 


. 8 8 


108 4 


. 8 7 


995 . 


5 6 


1017 


. 9 1 


1047 


. 1 


10 94 


. 5 2 


1131 


.7 


1132 


.41 


1129 


.72 


1141 


. 61 


1145 


. 6 


1121 


. 3 7 


107 6 


. 3 


1014 


. 62 


103 9 


. 7 5 


1071 


. 92 


10 97 


. 5 1 


1129 


.12 


1153 


.87 


1136 


. 98 


1135 


.12 


1126 


. 5 2 


1110 


. 8 9 


1082 


. 9 4 


1031 


. 3 2 


10 63 


. 3 3 


1100 


. 9 


1122 


.42 


1120 


. 62 


1141 


.09 


1154 


.8 


1151 


. 3 


1141 


. 4 3 


1125 


.42 


110 8 


. 2 8 


1049 


. 12 


1092 


. 55 


1124 


. 93 


1121 


. 58 


1123 


.57 


1154 


.33 


1193 


. 98 


1189 


. 47 


1174 


. 21 


1134 


. 22 


1103 


. 21 


1067 


. 64 


1113 


.28 


1130 


.49 


1118 


.21 


1125 


.38 


1160 


.22 


1234 


.75 


1240 


.27 


1200 


.85 


1141 


.13 


1104 


.25 


1076 


.79 


1109 


.31 


1117 


.12 


1114 


.58 


1123 


. 91 


1155 


.04 


1218 


.19 


1238 


.01 


1187 


.94 


1148 


.21 


1116 


.37 


1075 


. 66 


1092 


.29 


1113 


.24 


1114 


. 93 


1119 


.26 


1139 




1172 


.23 


1187 


.71 


1164 


.77 


1147 


.9 


1114 


. 95 


1053 


. 67 


1081 


.83 


1109 


.71 


1104 


.79 


1109 


. 94 


1120 


.85 


1141 


. 69 


1148 


. 69 


1128 


.75 


1103 


.89 


1073 


. 7 


1038 


. 97 


1069 


.31 


1102 


.85 


1097 


.81 


1100 


. 67 


1098 


.26 


1102 


.28 


1103 


.48 


1083 


.33 


1063 


.3 


1041 


. 13 



FIG. 6: a) no-error, b) 0.1% error added and c) 0.5% error added to G(r,r;cj). In this example, e(k) = —2(Cos[k x ] + Cos[k y ]) and 
rj(k) = 0.1 and uo — —1 in unit of t. 



Phase reconstruction algorithm 

Here, we write down the flowchart for the phase reconstruction algorithm that is followed to fix the phase of Go(R; w) which 
we get by taking the square root of the equation 

G (R; u)> = G(r - r; U) ~ G ^ = °J * G (R = 0; c)) 2 . (23) 

G(r = 0,r = 0;uj)-G (R = 0;uj) 

Upon taking the square root, we get Go(R; u) upto a phase of e Z7r . To fix this phase, we note that since the propagator in the 
continuum Go(R; uS) has to be a smooth well-behaved function for R ^ if it is to satisfy the Green's function equations of 
motion for the Hamiltonian operator and therefore its phase should also be a smooth and well-behaved as a function of R. To 
see this we start with the equation of motion for the non-interacting case in the continuum is 

{id I dt + V 2 /2m)G non (x, t; x' \ t') = S(x - x')S(t - t') (24) 

which upon Fourier transforming with respect to time gives 

(C + V 2 /2m)G non (x, x f ; Q = S(x - x f ) (25) 

For x 7^ x' , the above differential equation has no ill-behaved term and thus G non (x,x f ; Q has to be a well-behaved dif- 
ferentiate function. For the interacting case, the equations of motion is an infinite hierarchy of differential equations with 
the successive equations involving higher order Green's functions (See Vinay Ambegaokar's Chapter on The Green's Function 
Method in Superconductivity, Vol 1, edited by R. D. Parks). It is not clear to the author, how one could extend the non-interacting 
argument to this case. Instead, we argue as follows. As is usual in perturbation theory, the full propagator in momentum space 
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satisfies a Dyson's equation and is given by G(p, C) = (C — p 2 /2m — £(p; C)) _1 where E(p; C) is called the Self-energy and 
captures the effect of interactions. If this self-energy doesn't change the analytic structure of G(p, C) when compared to Gq(p, C) 
(More precisely, the pole at £ = p 2 /2m for the non-interaction survives, though it will get shifted off the real axis), then upon 
Fourier transforming to real space, the differentiability of G(x, x'\ Q will be preserved. Looking at Eq. [TTJs continuum version, 



V 2 G(R;uj) = Urn / dk— r — ^ (26) 



1 Z 4 ^ [fep^-R 

A7+ (2tt) 2 ./ BiZ< w + ^ _ \k\ 2 /2m - E(fe; cj) 

and if the pole structure of the integrand is same with and without E, then the differentiabiltity of Gq(x, x' £) implies differen- 
tiabiltity of G(x, a/; £). In the case of a lattice, Go(R; () is well-behaved for R ^ and at R = there is a kink in its phase 
(See the origin in Fig. |7]b) and d)). 

Similarly, Go(R;uj) 2, s phase should also be well-behaved as a function of space. This is condition that we impose on 
Go(R; uj) 2 while fixing phases. We start by making a spatial list of the phases as given by the Arg(z) function which restricts 
the phase obtained to the principal branch (— 7r, 7r]. Then, we start at the impurity site for which R = 0. As we move away from 
the origin, we multiply phase factors of e z2m7r to Go(R; uj) 2 = \Gq(R] a;) 2 |e^ priricipaZ for all R and the m's are so chosen that 
if \R'\ > \R\ then principa i + 2irm' > ^principal + 2irm. We implemented the choosing of m's in the following way : 

1) Define a monotoniser function that takes two arguments that lie between (— 7r,7r] and keeps adding 2tt to the second 
argument till it becomes greater than the first argument. mono(x, y) : Do y = y + 2tt Till y > x. 

The following steps are done in each of the symmetry-related octants in space and we write down the steps for the octant 
y > and x > y. 

2) Start at origin (0, 0). Move a step along x-axis to (x, y) = (1, 0). Then, mono((j) pr i nc i pa i(R =(x-l, y)), (^principal (R = 
(x,y)). 

3) Then do mono((j) pr i nc i P ai(R = (x, y)), (/^principal (R = (x, y + 1)) along y-direction till y = x. 

4) Move a step along x-axis. Do mono((j) pr i nc i pa i(R = {x — 1, y)), (^principal (R = (x, 0)) where the y of the first argument 
is highest integer such that x > (x — l) 2 + y 2 . 

5) Repeat step 3) and 4) till the whole octant is covered. 

Similar phase fixing is done for all the octants. Once this is done, the phase of Go(R; uj) is just half that of the phase-fixed 
Go(R; uj) 2 . Since, the phase of Go(R; uj) 2 has been made well-behaved, the phase of Gq(R; uj) will also be well-behaved which 
is what was desired. In Fig. U\ we show the result of doing the phase-fixing to numerically calculated Go(R;uj) 2 and also 
directly to Gq(R; uj) and find that they are in the correct ratio of two. 



Proof of self-energy relation 

In this section we prove that the relation between the electron and hole lifetimes, r]hoie(—uj) = — r] electron (uj). We will do 
this using the 2x2 Matsubara formalism. In this formalism, the Green's function for the non-lifetime broadened system in the 
normal state(i.e. no superconductivity) looks like 

**•>-' = C~ EW M, +e «°> < 27 > 

where uj n = (2n + 1)ttT is the fermionic Matsubara frequency. 

We will first prove the relation in the case of the normal electrons coupled to phonons. The self-energy in this case looks like 

Y,(k;iuj n ) = -— ^2 g(k - q, q)g(k, -q)D(q\ zO m )r 3 Go(A: - q; iuj n - zO m )r 3 

P,Qm 

(28) 

where Vt m = 2m7rT is the bosonic Matsubara frequency, T3 is the third componenet of Pauli matrices in the Nambu space, 
D(q; iQ m ) is the fourier-transform of the phonon's Green's function D(q; r) = — < T T [A(q; r)A(—q; 0)] > and evaluates to 

^^ ) = 5(sn^-in^W (29) 

where ft^ is the phonon dispersion. It has the following property : D(q; ift m ) = D(q; — iflm). The g(k,q) is the electron- 
phonon coupling strength coming from the electron-phonon interaction term 

Hel-ph = J^Y1 ^ C l + q,*%* A f- (30) 

k,q,a 
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FIG. 7: In these figures, we show the results of the Phase reconstruction algorithm. In a) and c), we show the phase (as a function of space 
around one quadrant of the impurity at origin) as evaluated using Arg(z) which restricts the values to one branch of the Argument function 
for Go(R; uj) and Go(R; uj) 2 respectively. In b) and d), we show the monotonised phase according to the algorithm described in this section. 
In e), we confirm that the ratio of the reconstructed phases of Go(R;uj) 2 and Go(R;uj) is identically two everywhere. In this example, 
e(k) = -2(Cos[k x ] + Cos[k y ]) and 77$) = 0.25 + 0.1(Cos[k x ] + Cos[k y \) andcj = -1. 
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Using the property D(<f; zO m ) = D(q; —ifL m ) and £l_ m = — f£ m , we can show that(suppressing momenta indices) 



E 22 (-i^ n ) ex ^ 



-icj n — + e — icj n + e — icj n — ifti 
-D(fi-i) -£>(fi ) -^(^l) 



iw n + — e iw n — e iw n + ifli — e 

-£>(fli) | -D(th) | -£>(»-i) 
iw n — ifli — e ioj n — e iui n — iCl-i — e 



— IL ill m c 

OC -Hii(^ n ) (31) 

Thus, by analytic continuation, £22 (— z) = — En(z) where H22 and En are the hole and electron self-energies respectively. 
Thus when we analytically continue till z — uj + iS where uj is real, we see that E 2 2(— uj — iS) — —En (a; + iS). From the 
analytic properties of Self-energy E(p; u±i5) — uj) =f cj)(see e.g., Eqn. 82 in Vinay Ambegaokar's Chapter on The 
Green's Function Method in Superconductivity, Vol 1, edited by R. D. Parks), we conclude that 

Vhole(-Uj) = -T]electron(uj) . . . QED (32) 

Also, the chemical potential shift is equal for both holes and electrons. This proof can be extended to higher orders in the 
electron-phonon coupling by noticing that all higher order terms contributing to self-energy contain odd number of fermion 
propagators, thus allowing the same kind of manipulation done above to go through analogously. This proof extends to other 
bosonic modes(e.g. spin wave modes) too since their propagators also satisfy D(q;iQ m ) = D(q; —iQ m ). This proof also 
extends to the case of lifetime broadening induced by electron-electron interaction by the same token that the self-energy terms 
always have odd number of fermion propagators. 
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